matlab
clc; clear; close all;

% 初始参数
g = 9.8;
v_missile = 300;
cloud_radius = 10;
cloud_duration = 20;
cloud_sink_speed = 3;

% 初始位置（以假目标为原点）
target_original = [0, 200, 0];
M1_start_original = [20000, 0, 2000];
FY1_start_original = [17800, 0, 1800];
FY2_start_original = [12000, 1400, 1400];
FY3_start_original = [6000, -3000, 700];

% 坐标转换
M1_start = M1_start_original - target_original;
FY1_start = FY1_start_original - target_original;
FY2_start = FY2_start_original - target_original;
FY3_start = FY3_start_original - target_original;

% 导弹速度方向
dir_missile = -M1_start / norm(M1_start);
v_missile_vec = v_missile * dir_missile;

% 定义目标函数（三架无人机各投放一枚烟幕弹）
objective_function = @(x) -calculate_multi_uav_obscuration_time(x, FY1_start, FY2_start, FY3_start, M1_start, v_missile_vec, g, cloud_radius, cloud_duration, cloud_sink_speed);

% 分层优化策略
fprintf('开始多无人机协同优化...\n');

% 第一阶段：单无人机单独优化
fprintf('第一阶段：单无人机单独优化\n');
single_opt_results = cell(3, 1);
uav_starts = {FY1_start, FY2_start, FY3_start};
uav_names = {'FY1', 'FY2', 'FY3'};

for i = 1:3
    fprintf('优化 %s...\n', uav_names{i});
    single_obj = @(x) -calculate_single_uav_obscuration_time(x, uav_starts{i}, M1_start, v_missile_vec, g, cloud_radius, cloud_duration, cloud_sink_speed);
    
    options_ga = optimoptions('ga', ...
        'PopulationSize', 50, ...
        'MaxGenerations', 30, ...
        'Display', 'off', ...
        'UseVectorized', false);
    
    lb = [0, 70, 0, 0];
    ub = [2*pi, 140, 20, 20];
    
    [x_opt, fval_opt] = ga(single_obj, 4, [], [], [], [], lb, ub, [], options_ga);
    single_opt_results{i} = x_opt;
    fprintf('%s 最优遮蔽时间: %.3f s\n', uav_names{i}, -fval_opt);
end

% 第二阶段：协同优化（以单机优化结果为初始点）
fprintf('\n第二阶段：多无人机协同优化\n');

% 构建初始解向量 [FY1_params, FY2_params, FY3_params]
x0 = [single_opt_results{1}, single_opt_results{2}, single_opt_results{3}];

% 设置优化选项
options_ps = optimoptions('patternsearch', ...
    'Display', 'iter', ...
    'FunctionTolerance', 1e-3, ...
    'StepTolerance', 1e-3, ...
    'MeshTolerance', 1e-3, ...
    'UseCompletePoll', true, ...
    'UseVectorized', false, ...
    'PlotFcn', {@psplotbestf, @psplotmeshsize});

% 定义变量边界
lb = repmat([0, 70, 0, 0], 1, 3);
ub = repmat([2*pi, 140, 20, 20], 1, 3);

% 运行模式搜索
[x_opt, fval_opt, exitflag_ps, output_ps] = patternsearch(objective_function, x0, [], [], [], [], lb, ub, [], options_ps);

% 提取最优参数
T_total_optimal = -fval_opt;

% 提取各无人机的参数
FY1_params = x_opt(1:4);
FY2_params = x_opt(5:8);
FY3_params = x_opt(9:12);

% 输出结果
fprintf('\n优化结果:\n');
fprintf('总遮蔽时间: %.3f s\n', T_total_optimal);

fprintf('\nFY1最优参数:\n');
fprintf('  飞行方向: %.2f 度\n', rad2deg(FY1_params(1)));
fprintf('  飞行速度: %.2f m/s\n', FY1_params(2));
fprintf('  投放时间: %.2f s\n', FY1_params(3));
fprintf('  起爆延迟: %.2f s\n', FY1_params(4));

fprintf('\nFY2最优参数:\n');
fprintf('  飞行方向: %.2f 度\n', rad2deg(FY2_params(1)));
fprintf('  飞行速度: %.2f m/s\n', FY2_params(2));
fprintf('  投放时间: %.2f s\n', FY2_params(3));
fprintf('  起爆延迟: %.2f s\n', FY2_params(4));

fprintf('\nFY3最优参数:\n');
fprintf('  飞行方向: %.2f 度\n', rad2deg(FY3_params(1)));
fprintf('  飞行速度: %.2f m/s\n', FY3_params(2));
fprintf('  投放时间: %.2f s\n', FY3_params(3));
fprintf('  起爆延迟: %.2f s\n', FY3_params(4));
% 计算并显示每架无人机的详细信息
[all_drop_points, all_burst_points, obscuration_times, total_obscuration] = calculate_multi_uav_details(...
    x_opt, FY1_start, FY2_start, FY3_start, M1_start, v_missile_vec, g, cloud_radius, cloud_duration, cloud_sink_speed);

for i = 1:3
    fprintf('\n%s:\n', uav_names{i});
    fprintf('  投放点坐标: (%.1f, %.1f, %.1f) m\n', all_drop_points(i, :));
    fprintf('  起爆点坐标: (%.1f, %.1f, %.1f) m\n', all_burst_points(i, :));
    fprintf('  单独遮蔽时间: %.3f s\n', obscuration_times(i));
end

fprintf('\n总遮蔽时间: %.3f s\n', total_obscuration);

% 可视化最优策略
visualize_multi_uav_solution(x_opt, FY1_start, FY2_start, FY3_start, M1_start, v_missile_vec, g, cloud_radius, cloud_duration, cloud_sink_speed);

% 保存结果到Excel文件
save_to_excel_result2(x_opt, all_drop_points, all_burst_points, obscuration_times, total_obscuration, uav_names);

% 单无人机遮蔽时间计算函数
function T = calculate_single_uav_obscuration_time(x, uav_start, M1_start, v_missile_vec, g, cloud_radius, cloud_duration, cloud_sink_speed)
    % 提取参数
    theta = x(1);
    v = x(2);
    t_drop = x(3);
    t_delay = x(4);
    t_burst = t_drop + t_delay;
    
    % 无人机速度向量
    v_fy = [v*cos(theta), v*sin(theta), 0];
    
    % 投放点位置
    drop_point = uav_start + v_fy * t_drop;
    
    % 起爆点位置
    g_vec = [0, 0, -g];
    burst_point = drop_point + v_fy * t_delay + 0.5 * g_vec * t_delay^2;
    
    % 导弹位置函数
    r_missile = @(t) M1_start + v_missile_vec * t;
    
    % 云团位置函数
    r_cloud = @(t) burst_point + [0, 0, -cloud_sink_speed] * (t - t_burst);
    
    % 确定时间范围
    t_start = t_burst;
    t_end = t_burst + cloud_duration;
    t_step = 0.01;
    t_vec = t_start:t_step:t_end;
    
    % 计算导弹和云团位置
    missile_pos = M1_start + v_missile_vec .* t_vec';
    
    cloud_pos = zeros(length(t_vec), 3);
    for i = 1:length(t_vec)
        t = t_vec(i);
        if t >= t_burst
            cloud_pos(i, :) = r_cloud(t);
        else
            cloud_pos(i, :) = [NaN, NaN, NaN];
        end
    end
    
    % 计算距离
    dists = sqrt(sum((missile_pos - cloud_pos).^2, 2));
    
    % 确定是否在有效范围内
    in_cloud = dists <= cloud_radius & t_vec' >= t_burst & t_vec' <= t_burst + cloud_duration;
    
    % 计算有效遮蔽时间
    T = sum(in_cloud) * t_step;
end

% 多无人机总遮蔽时间计算函数
function T_total = calculate_multi_uav_obscuration_time(x, FY1_start, FY2_start, FY3_start, M1_start, v_missile_vec, g, cloud_radius, cloud_duration, cloud_sink_speed)
    % 提取参数
    FY1_params = x(1:4);
    FY2_params = x(5:8);
    FY3_params = x(9:12);
    
    % 计算各无人机的烟幕弹参数
    uav_starts = {FY1_start, FY2_start, FY3_start};
    all_params = {FY1_params, FY2_params, FY3_params};
    
    burst_points = zeros(3, 3);
    t_bursts = zeros(3, 1);
    
    for i = 1:3
        theta = all_params{i}(1);
        v = all_params{i}(2);
        t_drop = all_params{i}(3);
        t_delay = all_params{i}(4);
        t_burst = t_drop + t_delay;
        
        % 无人机速度向量
        v_fy = [v*cos(theta), v*sin(theta), 0];
        
        % 投放点位置
        drop_point = uav_starts{i} + v_fy * t_drop;
        
        % 起爆点位置
        g_vec = [0, 0, -g];
        burst_points(i, :) = drop_point + v_fy * t_delay + 0.5 * g_vec * t_delay^2;
        t_bursts(i) = t_burst;
    end
    
     % 导弹位置函数
    r_missile = @(t) M1_start + v_missile_vec * t;
    
    % 云团位置函数
    r_cloud1 = @(t) burst_points(1, :) + [0, 0, -cloud_sink_speed] * (t - t_bursts(1));
    r_cloud2 = @(t) burst_points(2, :) + [0, 0, -cloud_sink_speed] * (t - t_bursts(2));
    r_cloud3 = @(t) burst_points(3, :) + [0, 0, -cloud_sink_speed] * (t - t_bursts(3));
    
    % 确定时间范围
    t_start = min(t_bursts);
    t_end = max(t_bursts) + cloud_duration;
    t_step = 0.01;
    t_vec = t_start:t_step:t_end;
    
    % 计算导弹和云团位置
    missile_pos = M1_start + v_missile_vec .* t_vec';
    
    cloud1_pos = zeros(length(t_vec), 3);
    cloud2_pos = zeros(length(t_vec), 3);
    cloud3_pos = zeros(length(t_vec), 3);
    
    for i = 1:length(t_vec)
        t = t_vec(i);
        if t >= t_bursts(1)
            cloud1_pos(i, :) = r_cloud1(t);
        else
            cloud1_pos(i, :) = [NaN, NaN, NaN];
        end
        
        if t >= t_bursts(2)
            cloud2_pos(i, :) = r_cloud2(t);
        else
            cloud2_pos(i, :) = [NaN, NaN, NaN];
        end
        
        if t >= t_bursts(3)
            cloud3_pos(i, :) = r_cloud3(t);
        else
            cloud3_pos(i, :) = [NaN, NaN, NaN];
        end
    end
    
    % 计算距离
    dist1 = sqrt(sum((missile_pos - cloud1_pos).^2, 2));
    dist2 = sqrt(sum((missile_pos - cloud2_pos).^2, 2));
    dist3 = sqrt(sum((missile_pos - cloud3_pos).^2, 2));
    
    % 确定是否在有效范围内
    in_cloud1 = dist1 <= cloud_radius & t_vec' >= t_bursts(1) & t_vec' <= t_bursts(1) + cloud_duration;
    in_cloud2 = dist2 <= cloud_radius & t_vec' >= t_bursts(2) & t_vec' <= t_bursts(2) + cloud_duration;
    in_cloud3 = dist3 <= cloud_radius & t_vec' >= t_bursts(3) & t_vec' <= t_bursts(3) + cloud_duration;
    
    % 计算总遮蔽时间（并集）
    in_any_cloud = in_cloud1 | in_cloud2 | in_cloud3;
    T_total = sum(in_any_cloud) * t_step;
end

% 计算多无人机详细信息
function [all_drop_points, all_burst_points, obscuration_times, total_obscuration] = calculate_multi_uav_details(...
    x, FY1_start, FY2_start, FY3_start, M1_start, v_missile_vec, g, cloud_radius, cloud_duration, cloud_sink_speed)
    
    % 提取参数
    FY1_params = x(1:4);
    FY2_params = x(5:8);
    FY3_params = x(9:12);
    
    % 计算各无人机的烟幕弹参数
    uav_starts = {FY1_start, FY2_start, FY3_start};
    all_params = {FY1_params, FY2_params, FY3_params};
    
    all_drop_points = zeros(3, 3);
    all_burst_points = zeros(3, 3);
    t_bursts = zeros(3, 1);
    
    for i = 1:3
        theta = all_params{i}(1);
        v = all_params{i}(2);
        t_drop = all_params{i}(3);
        t_delay = all_params{i}(4);
        t_burst = t_drop + t_delay;
        
        % 无人机速度向量
        v_fy = [v*cos(theta), v*sin(theta), 0];
        
        % 投放点位置
        all_drop_points(i, :) = uav_starts{i} + v_fy * t_drop;
        
        % 起爆点位置
        g_vec = [0, 0, -g];
        all_burst_points(i, :) = all_drop_points(i, :) + v_fy * t_delay + 0.5 * g_vec * t_delay^2;
        t_bursts(i) = t_burst;
    end
    
    % 导弹位置函数
    r_missile = @(t) M1_start + v_missile_vec * t;
    
    % 云团位置函数
    r_cloud1 = @(t) all_burst_points(1, :) + [0, 0, -cloud_sink_speed] * (t - t_bursts(1));
    r_cloud2 = @(t) all_burst_points(2, :) + [0, 0, -cloud_sink_speed] * (t - t_bursts(2));
    r_cloud3 = @(t) all_burst_points(3, :) + [0, 0, -cloud_sink_speed] * (t - t_bursts(3));
    
    % 计算每枚烟幕弹的单独遮蔽时间
    obscuration_times = zeros(3, 1);
    t_step = 0.001;
    
    for i = 1:3
        t_burst = t_bursts(i);
        t_start = t_burst;
        t_end = t_burst + cloud_duration;
        t_vec = t_start:t_step:t_end;
        
        if i == 1
            r_cloud = r_cloud1;
        elseif i == 2
            r_cloud = r_cloud2;
        else
            r_cloud = r_cloud3;
        end
        
        dists = arrayfun(@(t) norm(r_missile(t) - r_cloud(t)), t_vec);
        in_cloud = dists <= cloud_radius;
        obscuration_times(i) = sum(in_cloud) * t_step;
    end
    
    % 计算总遮蔽时间
    t_start = min(t_bursts);
    t_end = max(t_bursts) + cloud_duration;
    t_vec = t_start:t_step:t_end;
    
    in_any_cloud = false(length(t_vec), 1);
    for i = 1:length(t_vec)
        t = t_vec(i);
        if t >= t_bursts(1) && t <= t_bursts(1) + cloud_duration
            in_any_cloud(i) = in_any_cloud(i) || (norm(r_missile(t) - r_cloud1(t)) <= cloud_radius);
        end
        if t >= t_bursts(2) && t <= t_bursts(2) + cloud_duration
            in_any_cloud(i) = in_any_cloud(i) || (norm(r_missile(t) - r_cloud2(t)) <= cloud_radius);
        end
        if t >= t_bursts(3) && t <= t_bursts(3) + cloud_duration
            in_any_cloud(i) = in_any_cloud(i) || (norm(r_missile(t) - r_cloud3(t)) <= cloud_radius);
        end
    end
    
    total_obscuration = sum(in_any_cloud) * t_step;
end
% 可视化多无人机解决方案
function visualize_multi_uav_solution(x, FY1_start, FY2_start, FY3_start, M1_start, v_missile_vec, g, cloud_radius, cloud_duration, cloud_sink_speed)
    % 提取参数
    FY1_params = x(1:4);
    FY2_params = x(5:8);
    FY3_params = x(9:12);
    
    % 计算各无人机的烟幕弹参数
    uav_starts = {FY1_start, FY2_start, FY3_start};
    all_params = {FY1_params, FY2_params, FY3_params};
    uav_names = {'FY1', 'FY2', 'FY3'};
    colors = {'r', 'g', 'b'};
    
    all_drop_points = zeros(3, 3);
    all_burst_points = zeros(3, 3);
    t_bursts = zeros(3, 1);
    
    for i = 1:3
        theta = all_params{i}(1);
        v = all_params{i}(2);
        t_drop = all_params{i}(3);
        t_delay = all_params{i}(4);
        t_burst = t_drop + t_delay;
        
        % 无人机速度向量
        v_fy = [v*cos(theta), v*sin(theta), 0];
        
        % 投放点位置
        all_drop_points(i, :) = uav_starts{i} + v_fy * t_drop;
        
        % 起爆点位置
        g_vec = [0, 0, -g];
        all_burst_points(i, :) = all_drop_points(i, :) + v_fy * t_delay + 0.5 * g_vec * t_delay^2;
        t_bursts(i) = t_burst;
    end
    
    % 导弹位置函数
    r_missile = @(t) M1_start + v_missile_vec * t;
    
    % 云团位置函数
    r_cloud1 = @(t) all_burst_points(1, :) + [0, 0, -cloud_sink_speed] * (t - t_bursts(1));
    r_cloud2 = @(t) all_burst_points(2, :) + [0, 0, -cloud_sink_speed] * (t - t_bursts(2));
    r_cloud3 = @(t) all_burst_points(3, :) + [0, 0, -cloud_sink_speed] * (t - t_bursts(3));
    
    % 确定时间范围
    t_start = min(t_bursts) - 5;
    t_end = max(t_bursts) + cloud_duration + 5;
    t_step = 0.01;
    t_vec = t_start:t_step:t_end;
    
    % 计算导弹和云团位置
    missile_pos = zeros(length(t_vec), 3);
    cloud1_pos = zeros(length(t_vec), 3);
    cloud2_pos = zeros(length(t_vec), 3);
    cloud3_pos = zeros(length(t_vec), 3);
    
    for i = 1:length(t_vec)
        t = t_vec(i);
        missile_pos(i, :) = r_missile(t);
        
        if t >= t_bursts(1)
            cloud1_pos(i, :) = r_cloud1(t);
        else
            cloud1_pos(i, :) = [NaN, NaN, NaN];
        end
        
        if t >= t_bursts(2)
            cloud2_pos(i, :) = r_cloud2(t);
        else
            cloud2_pos(i, :) = [NaN, NaN, NaN];
        end
        
        if t >= t_bursts(3)
            cloud3_pos(i, :) = r_cloud3(t);
        else
            cloud3_pos(i, :) = [NaN, NaN, NaN];
        end
    end
    
    % 计算距离
    dist1 = sqrt(sum((missile_pos - cloud1_pos).^2, 2));
    dist2 = sqrt(sum((missile_pos - cloud2_pos).^2, 2));
    dist3 = sqrt(sum((missile_pos - cloud3_pos).^2, 2));
    
    % 确定是否在有效范围内
    in_cloud1 = dist1 <= cloud_radius & t_vec' >= t_bursts(1) & t_vec' <= t_bursts(1) + cloud_duration;
    in_cloud2 = dist2 <= cloud_radius & t_vec' >= t_bursts(2) & t_vec' <= t_bursts(2) + cloud_duration;
    in_cloud3 = dist3 <= cloud_radius & t_vec' >= t_bursts(3) & t_vec' <= t_bursts(3) + cloud_duration;
    
    in_any_cloud = in_cloud1 | in_cloud2 | in_cloud3;
    
    % 创建可视化
    figure('Position', [100, 100, 1400, 900], 'Name', '多无人机烟幕干扰策略分析');
    
    % 子图1: 距离随时间变化
    subplot(2,2,1);
    plot(t_vec, dist1, 'r-', 'LineWidth', 1);
    hold on;
    plot(t_vec, dist2, 'g-', 'LineWidth', 1);
    plot(t_vec, dist3, 'b-', 'LineWidth', 1);
    plot(t_vec, cloud_radius*ones(size(t_vec)), 'k--', 'LineWidth', 1.5);
    
    % 标记有效遮蔽区域
    area_x = t_vec(in_any_cloud);
    if ~isempty(area_x)
        area_y = zeros(size(area_x)) + max([max(dist1), max(dist2), max(dist3)]);
        fill([min(area_x), min(area_x), max(area_x), max(area_x)], [0, max(area_y), max(area_y), 0], 'y', 'FaceAlpha', 0.2, 'EdgeColor', 'none');
    end
    
    legend('FY1烟幕距离', 'FY2烟幕距离', 'FY3烟幕距离', '有效半径阈值', '有效遮蔽时段', 'Location', 'best');
    xlabel('时间 (s)');
    ylabel('导弹与云团中心距离 (m)');
    title('导弹与各烟幕云团中心距离随时间变化');
    grid on;
    
    % 子图2: 是否在有效范围内
    subplot(2,2,2);
    stairs(t_vec, in_cloud1, 'r-', 'LineWidth', 1.5);
    hold on;
    stairs(t_vec, in_cloud2, 'g-', 'LineWidth', 1.5);
    stairs(t_vec, in_cloud3, 'b-', 'LineWidth', 1.5);
    stairs(t_vec, in_any_cloud, 'k-', 'LineWidth', 2);
    
    legend('FY1', 'FY2', 'FY3', '总有效遮蔽', 'Location', 'best');
    xlabel('时间 (s)');
    ylabel('是否在有效范围内');
    title('导弹是否在各烟幕云团有效范围内');
    ylim([-0.1, 1.1]);
    grid on;
    
    % 子图3: 三维轨迹可视化
    subplot(2,2,3);
    % 绘制导弹轨迹
    plot3(missile_pos(:,1), missile_pos(:,2), missile_pos(:,3), 'k-', 'LineWidth', 2);
    hold on;
    
    % 绘制无人机初始位置
    plot3(FY1_start(1), FY1_start(2), FY1_start(3), 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r');
    plot3(FY2_start(1), FY2_start(2), FY2_start(3), 'go', 'MarkerSize', 8, 'MarkerFaceColor', 'g');
    plot3(FY3_start(1), FY3_start(2), FY3_start(3), 'bo', 'MarkerSize', 8, 'MarkerFaceColor', 'b');
    
    % 绘制云团中心轨迹
    valid_cloud1 = ~isnan(cloud1_pos(:,1));
    plot3(cloud1_pos(valid_cloud1,1), cloud1_pos(valid_cloud1,2), cloud1_pos(valid_cloud1,3), 'r--', 'LineWidth', 1.5);
    
    valid_cloud2 = ~isnan(cloud2_pos(:,1));
    plot3(cloud2_pos(valid_cloud2,1), cloud2_pos(valid_cloud2,2), cloud2_pos(valid_cloud2,3), 'g--', 'LineWidth', 1.5);
    
    valid_cloud3 = ~isnan(cloud3_pos(:,1));
    plot3(cloud3_pos(valid_cloud3,1), cloud3_pos(valid_cloud3,2), cloud3_pos(valid_cloud3,3), 'b--', 'LineWidth', 1.5);
    
    % 标记关键点
    for i = 1:3
        plot3(all_drop_points(i,1), all_drop_points(i,2), all_drop_points(i,3), [colors{i} 'o'], 'MarkerSize', 8, 'MarkerFaceColor', colors{i});
        plot3(all_burst_points(i,1), all_burst_points(i,2), all_burst_points(i,3), [colors{i} 'x'], 'MarkerSize', 8, 'LineWidth', 2);
    end
    
    plot3(0, 0, 0, 'kx', 'MarkerSize', 10, 'LineWidth', 2); % 假目标位置
    
    legend('导弹轨迹', 'FY1初始位置', 'FY2初始位置', 'FY3初始位置', ...
           'FY1云团', 'FY2云团', 'FY3云团', ...
           'FY1投放点', 'FY1起爆点', 'FY2投放点', 'FY2起爆点', 'FY3投放点', 'FY3起爆点', '假目标', 'Location', 'best');
    xlabel('X (m)');
    ylabel('Y (m)');
    zlabel('Z (m)');
    title('三维轨迹可视化');
    grid on;
    view(3);
    
    % 子图4: 烟幕弹时序图
    subplot(2,2,4);
    % 绘制烟幕弹投放和起爆时间线
    y_values = [3, 2, 1]; % FY1在最上面
    
    for i = 1:3
        t_drop = all_params{i}(3);
        t_delay = all_params{i}(4);
        t_burst = t_drop + t_delay;
        
        % 投放时间
        plot([t_drop, t_drop], [y_values(i)-0.4, y_values(i)+0.4], 'k-', 'LineWidth', 2);
        hold on;
        text(t_drop, y_values(i)+0.5, sprintf('%s投放', uav_names{i}), 'HorizontalAlignment', 'center');
        
        % 起爆时间
        plot([t_burst, t_burst], [y_values(i)-0.4, y_values(i)+0.4], [colors{i} '-'], 'LineWidth', 2);
        text(t_burst, y_values(i)+0.5, sprintf('%s起爆', uav_names{i}), 'HorizontalAlignment', 'center');
        
        % 烟幕持续时间
        plot([t_burst, t_burst+cloud_duration], [y_values(i), y_values(i)], [colors{i} '-'], 'LineWidth', 3);
    end
    
    xlabel('时间 (s)');
    ylabel('无人机');
    ylim([0, 4]);
    set(gca, 'YTick', [1, 2, 3], 'YTickLabel', {'FY3', 'FY2', 'FY1'});
    title('各无人机烟幕弹投放和起爆时序');
    grid on;
end

% 保存结果到Excel文件（result2.xlsx）
function save_to_excel_result2(x_opt, all_drop_points, all_burst_points, obscuration_times, total_obscuration, ~)
    % 创建数据表
    data = cell(4, 10);
    
    % 表头
    data{1, 1} = '无人机编号';
    data{1, 2} = '无人机运动方向';
    data{1, 3} = '无人机运动速度 (m/s)';
    data{1, 4} = '烟幕干扰弹投放点的x坐标 (m)';
    data{1, 5} = '烟幕干扰弹投放点的y坐标 (m)';
    data{1, 6} = '烟幕干扰弹投放点的z坐标 (m)';
    data{1, 7} = '烟幕干扰弹起爆点的x坐标 (m)';
    data{1, 8} = '烟幕干扰弹起爆点的y坐标 (m)';
    data{1, 9} = '烟幕干扰弹起爆点的z坐标 (m)';
    data{1, 10} = '有效干扰时长 (s)';
    
    % 提取参数
    FY1_params = x_opt(1:4);
    FY2_params = x_opt(5:8);
    FY3_params = x_opt(9:12);
    
    % FY1数据
    data{2, 1} = 'FY1';
    data{2, 2} = rad2deg(FY1_params(1)); % 转换为角度
    data{2, 3} = FY1_params(2);
    data{2, 4} = all_drop_points(1, 1);
    data{2, 5} = all_drop_points(1, 2);
    data{2, 6} = all_drop_points(1, 3);
    data{2, 7} = all_burst_points(1, 1);
    data{2, 8} = all_burst_points(1, 2);
    data{2, 9} = all_burst_points(1, 3);
    data{2, 10} = obscuration_times(1);
    
    % FY2数据
    data{3, 1} = 'FY2';
    data{3, 2} = rad2deg(FY2_params(1));
    data{3, 3} = FY2_params(2);
    data{3, 4} = all_drop_points(2, 1);
    data{3, 5} = all_drop_points(2, 2);
    data{3, 6} = all_drop_points(2, 3);
    data{3, 7} = all_burst_points(2, 1);
    data{3, 8} = all_burst_points(2, 2);
    data{3, 9} = all_burst_points(2, 3);
    data{3, 10} = obscuration_times(2);
    
    % FY3数据
    data{4, 1} = 'FY3';
    data{4, 2} = rad2deg(FY3_params(1));
    data{4, 3} = FY3_params(2);
    data{4, 4} = all_drop_points(3, 1);
    data{4, 5} = all_drop_points(3, 2);
    data{4, 6} = all_drop_points(3, 3);
    data{4, 7} = all_burst_points(3, 1);
    data{4, 8} = all_burst_points(3, 2);
    data{4, 9} = all_burst_points(3, 3);
    data{4, 10} = obscuration_times(3);
    
    % 写入Excel文件
    filename = 'result2.xlsx';
    writecell(data, filename, 'Sheet', 'Sheet1');
    
    fprintf('结果已保存到 %s\n', filename);
    
    % 添加备注
   e = {
        '注：以x轴为正向，逆时针方向为正，取值0~360（度）。';
        sprintf('总有效遮蔽时间: %.3f s', total_obscuration);
    };
    
    % 尝试将备注添加到Excel文件
    try
        excel = actxserver('Excel.Application');
        excel.Visible = 0;
        workbook = excel.Workbooks.Open([pwd, '\', filename]);
        sheets = workbook.Sheets;
        sheet = sheets.Item('Sheet1');
       e;
        % 在适当位置添加备注
        sheet.Range('A6').Value = e{1};
        sheet.Range('A7').Value = e{2};
        
        workbook.Save;
        workbook.Close;
        excel.Quit;
    catch
        fprintf('无法添加备注到Excel文件，请手动添加:\n');
        fprintf('%s\n', e{1});
        fprintf('%s\n', e{2});
    end
end